rm(list = ls())
# load all required libraries
# Analysis was conducted in R version 4.2.1
library(dplyr)
library(sjPlot)
library(tidyr)
library(mice)
library(psych)
library(MASS)
library(stargazer)
library(ordinal)
library(openxlsx)

############party identification
data_tw$party <- data_tw$pid
data_tw$party[data_tw$pid %in% c(1,4,6)] <- "Pan Blue"
data_tw$party[data_tw$pid %in% c(2,5,7,8,9)] <- "Pan Green"
data_tw$party[data_tw$pid == 3] <- "TPP"
data_tw$party[data_tw$pid %in% 10:11] <- "Non Partisans"
data_tw$party <- as.factor(data_tw$party)
table(data_tw$party)
100*prop.table(table(data_tw$party))
partyTable<-table(data_tw$party)
barplot(partyTable, 
        main="Barplot of Party Identification",
        xlab="Party",
        col = "blue")
###############gender
data_tw$genderR<-factor(data_tw$gender, 
                             levels=c(1,2), 
                             labels=c("female", "Male"))
sexTable<-table(data_tw$genderR)
prop.table(table(data_tw$genderR))
###49,53% are female
###50,46% male
barplot(sexTable, 
        main="Barplot of sex",
        xlab="Sex",
        col = "yellow")
#Histogram of age with formatting
table(data_tw$age)
data_tw$age[data_tw$age=="4"]<-18 ###to remove outliner
data_tw$age[data_tw$age=="15"]<-18
data_tw$age[data_tw$age=="17"]<-18
data_tw$age[data_tw$age>=18&data_tw$age<30]<-1
data_tw$age[data_tw$age>=30&data_tw$age<50]<-2
data_tw$age[data_tw$age>=50&data_tw$age<70]<-3
data_tw$age<-as.factor(data_tw$age)
hist(data_tw$age,
     breaks=10,
     xlab="Age in years", 
     main="Histogram of Age", 
     col="lightblue")
summary(data_tw$age)
####### Minimum age is 18.
##Maximum age is 69
##25% of cases are under 29 years of age
###50% of cases are under 37 years of age
###75% of cases are under 45 years of age
######description
data_tw_desc<-subset(data_tw,
                      select=(c("ai_priv",
                                "Combined_reg",
                                "ri_military",
                                "ai_risks",
                                "pid",
                                "age",
                                "gender",
                                "edu")))
ai_tw_description<-psych::describe(data_tw_desc)
write.xlsx(ai_tw_description,"ai_tw_description.xlsx")
#########education
data_tw$edu_binary <- ifelse(data_tw$edu != 5, "0 - no Master", "1 - Master or higher")
data_tw$edu_binary <- as.factor(data_tw$edu_binary)
table(data_tw$edu_binary)
data_tw$eduR<-ordered(data_tw$edu,
                      levels=c(1,2,3,4,5), 
                      labels=c("Elementary School", 
                               "Junior High School",
                               "High School",
                               "University",
                               "Master or Higher"))
100*prop.table(table(data_tw$eduR))
EducTable <- table(data_tw$eduR)
par(las=2)
barplot(EducTable,
        main="Barplot of Education", 
        col="red",
        cex.names=0.5)
#side-by-side boxplot (Gender ~ Age)
boxplot(data_tw$age ~ data_tw$genderR,
        main="Boxplot of Age by Sex",
        ylab="Age in years",
        xlab="Sex",
        col="darkgreen",
        notch=TRUE)

#Creating a dataframe that includes the relevant items only 
data_tw_alpha<-subset(data_tw,
                           select=(c("ai_priv",
                                     "reg_ban",
                                     "pid",
                                     "age",
                                     "gender")))


#Missing data plot using the visdat package
visdat::vis_miss(data_tw_alpha)


#######
table(data_tw$ai_priv)
data_tw$ai_priv<-ordered(data_tw$ai_priv, levels = 1:5)
data_tw$reg_ban<-ordered(data_tw$reg_ban, levels = 1:7)
data_tw$Combined_reg<-ordered(data_tw$Combined_reg, levels = 1:7)
data_tw$pid <- as.numeric(data_tw$pid)
data_tw$edu <- as.numeric(data_tw$edu)
data_tw$age<- as.numeric(data_tw$age)
data_tw$ai_prof<- as.numeric(data_tw$ai_prof)
data_tw$ai_priv<- as.numeric(data_tw$ai_priv)
data_tw$Combined_reg<-as.numeric(data_tw$Combined_reg)
data_tw$ri_military<-as.numeric(data_tw$ri_military)
data_tw$ai_risks<-as.numeric(data_tw$ai_risks)
data_tw$polor_con<-as.numeric(data_tw$polor_con)
data_tw$reg_ban<-as.numeric(data_tw$reg_ban)

library(ordinal)
data_tw$ai_risks <- rowMeans(data_tw[ ,c("ri_job", "ri_military", "ri_control")])
data_tw$Combined_reg <- rowMeans(data_tw[, c("reg_state", "reg_plat","reg_ban")], na.rm = TRUE)
data_tw$Combined_AI_Use <- rowMeans(data_tw[, c("ai_prof", "ai_priv")], na.rm = TRUE)
m_new<-clm(Combined_reg ~ ai_priv + pid + age + gender + edu, data = data_tw)
m_new_mil<-clm(Combined_reg ~ ai_risks + pid + age + gender + edu, data = data_tw)
m_general <- clm(reg_ban ~ ai_priv + age + gender + edu, data = data_tw) ###combines private and professional usage of AI
m_1 <- clm(reg_ban ~ pid + age + gender + edu, data = data_tw) ###combines private and professional usage of AI
m_2 <- clm(reg_ban ~ ai_prof + age + gender + edu, data = data_tw) 
m_3 <- clm(reg_ban ~ pid + age + gender + edu, data = data_tw) 
summary(m_new)
summary(m_new_mil)
summary(m4)
library(stargazer)
stargazer(m4, type="html", ord.intercepts = TRUE,out="taiwan_ai.html")
rm(list = ls())
ctable <- coef(summary(m4))
t_vals <- ctable[, grep("t", colnames(ctable), ignore.case = TRUE)]
p_values <- pnorm(abs(t_vals), lower.tail = FALSE) * 2
ctable <- cbind(ctable, "p value" = p_values)
print(ctable)

##########deneme
data_tw$ai_priv<-ordered(data_tw$ai_priv, levels = 1:5)
m_deneme<-clm(ai_priv ~ Combined_reg + party + age + gender + edu, data = data_tw)
m_deneme_1<-clm(ai_priv ~ ri_military + party + age + gender + edu, data = data_tw)

summary(m_deneme)
summary(m_deneme_1)

library(stargazer)
stargazer(m_deneme,m_deneme_1, type="html", ord.intercepts = TRUE,out="taiwan_ai_deneme.html")

################FINAL MODEL
library(ordinal)
data_tw$Combined_reg <- rowMeans(data_tw[, c("reg_state", "reg_plat","reg_ban")], na.rm = TRUE)
data_tw$reg_ban<-ordered(data_tw$reg_ban, levels = 1:7)
data_tw$Combined_reg<-ordered(data_tw$Combined_reg, levels = 1:7)
data_tw$pid <- as.numeric(data_tw$pid)
data_tw$edu <- as.numeric(data_tw$edu)
data_tw$age<- as.numeric(data_tw$age)
data_tw$ai_priv<- as.numeric(data_tw$ai_priv)
m_new<-clm(Combined_reg ~ ai_priv + pid + age + gender + edu, data = data_tw)
m_new_1<-clm(Combined_reg ~ ri_military + pid + age + gender + edu, data = data_tw)
summary(m_new)
##############compare models
AIC(m_1, m_2, m_general)
# Perform Likelihood Ratio Test
lrtest_result <- anova(m_1, m_2, m_general)
print(lrtest_result)


#############Check Model Fit
install.packages("pscl")
library(pscl)
pR2(m_deneme)
####McFadden's pseudo r-squared result is 0.18
AIC(m_deneme)
###AIC result is low
pR2(m_deneme_1)
#######McFadden's pseudo r-squared result is 0.19
AIC(m_deneme_1)
###AIC result is low

########Multicollinearity
library(car)
vif(m_deneme)

